Vacuum Polarization and Screening of Supercritical Impurities in Graphene 
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Screening of charge impurities in graphene is analyzed using the exact solution for vacuum po- 
larization obtained from the massless Dirac-Kepler problem. For the impurity charge below certain 
critical value no density perturbation is found away from the impurity, in agreement with the linear 
response theory result. For supercritical charge, however, the polarization distribution is shown to 
have a power law profile, leading to screening of the excess charge at large distances. The Dirac- 
Kepler scattering states give rise to standing wave oscillations in the local density of states which 
appear and become prominent in the supercritical regime. 



Massless Dirac excitations in graphene [l| provide 
an interesting realization of quantum electrodynamics 
(QED) in dimension two [2]. Because of zero mass and 
strong interactions, characterized by a large "fine struc- 
ture constant" a = /hvp ^ 2.5, where vp ^ 10^ m/s 
is the Fermi velocity, this material breaks away from 
the perturbative QED paradigm. One of the phenom- 
ena fundamental in QED, expected to become ultra- 
strong in graphene, is "vacuum polarization" induced by 
charge impurities. Scattering on charge impurities fea- 
tures prominently in transport properties of graphene, 
where it is believed to be the leading factor limiting elec- 
tron mobilityjpl, y, l5|] , providing an explanation for the 
conductivity [61] linear in the carrier concentration. Al- 
though the problem of Coulomb scattering by charge im- 
purities received a lot of attention [I, [J, |5|, III, l8|, , the 
key question of screening of the impurity potential out- 
side the weak coupling regime has not been adequately 
addressed [1, 0] . 

Here we present an accurate nonperturbative treat- 
ment of this problem based on the vacuum polariza- 
tion found from the exact solution of the 2d Dirac- 
Kepler problem. There are two qualitatively different 
regimes emerging from this solution, which are somewhat 
analogous to those known in QED of heavy and super- 
heavy atoms jTlj . The Coulomb potential of subcritical 
strength can be treated as a mathematical singularity in 
solving the Dirac equation, while in the supercritical case 
a consistent solution is only possible after finite radius of 
charge distribution in a nucleus is accounted for [13] • We 
shall see that a similar phenomenon takes place in our 
problem at the critical charge value 
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where hi is the effective dielectric constant. For the case 
when screening is solely due to the graphene electrons, 
the RPA approach [2] gives /^rpa ^ 5. With e^/hvp ^ 
2.5, this yields a critical value Zc ^ 1. 

The most prominent effect in our problem, arising at 
supercritical is the change in the character of polar- 



ization of the Dirac vacuum. While at < ^ the polar- 
ization charge Qp is localized on the scale of the impurity 
radius, exhibiting no long range tail [9], for supercritical 
/3 the solution of the massless Dirac equation predicts 
a power law for the spatial profile of polarization. For 
^ < /3 < |, when just the lowest angular momentum 
channels of the Dirac equation are overcritical, we find 
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where A/" = 4 is the combined spin and valley degener- 
acy of graphene. The result ([2]), valid for noninter acting 
fermions, is somewhat modified at higher [3 (see Eq. (|29|) ). 

The result ([2]) can be used to describe screening in 
an interacting system in a selfconsistent renormalization 
group (RG) fashion. The KG flow for polarization cloud 
is constructed by proceeding from the lattice scale p = vq 
to larger p, treating the net polarization charge within 
radius p as an effective point charge and using it to 
determine polarization at larger distances. As a result, 
the net charge P{p) flows from its initial value P{ro) to 
lower values at larger distances. The net polarization 
charge ([2|) within the annulus pi < p < p2 equals 5Z = 
—A" sign /^^ ln(p2/pi), which leads to the RG equation 
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Integrating the flow ([3|) we find that it terminates at a 
distance p^ = tq exp ( ^'^7^^ cosh" "'^(2/3)) where P reaches 
the critical value ([1]). In contrast to screening in metals, 
here the polarization build-up brings the net charge down 
to the critical value Pc that remains unscreened at larger 
distances p ^ p^. The RG treatment is applicable when 
the RG flow is slow, i.e., when the right-hand side of 
Eq. ([3|) is small. Thus the RG framework is adequate 
near the criticality, p ^ pc^ where 7 is small, even in the 
case of strong coupling, e'^/tvhvF ~ 1, and the predicted 
termination of screening at large p is universal. 

Our treatment of vacuum polarization relies on the ex- 
act solution of the Dirac-Kepler problem from which we 
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extract scattering phases and use them in the Friedel 
sum rule framework to evaluate the screening charge. 
The phases are found to behave differently for /3 < /3c 
and P > Pc In the first case, the essential part 
of the phase is p\nkp^ while in the second case it is 
p\nkp — 7 sign In Zero- It is important to realize that 
the term plnkp is the same for all angular momentum 
channels. Such a contribution to the phase, as Ref.[12| 
insightfully remarks, arises from quasiclassical dynam- 
ics at large distances, and thus has nothing to do with 
scattering. In agreement with this, we find that it does 
not contribute to polarization at finite p, while the term 
—7 sign/? In Zero gives rise to the power law in ([2]). 

Now we turn to the analysis of the massless Dirac equa- 
tion in dimension two in a central potential V{p): 
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Introducing polar coordinates x-\-iy = pe*'^, we separate 
angular harmonics of the two-component wave function 
and seek the solution in the form [2^] 
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with integer angular quantum number m. The terms w 
and V represent the incoming and outgoing waves. The 
parameters s and k are determined by the behavior at 
small and large p. For V{p) = Ze^ j p^ we find 



{ra^\f-ii\ k = -j^, (6) 

where P is the dimensionless coupling ([T]). (The minus 
sign in (|6|) is chosen to make /c > in the Fermi sea.) 
The solution ([5]) behaves differently for \p\ < |m + ||, 
when the exponent s is real, and \P\ > |m + ^|, when s 
becomes complex imaginary. 

The ansatz ([5]), substituted into Dirac equation, yields 
coupled equations for the functions w{p) and v{p): 

pdpW + (5 + + 2ikp) w - {m^\)v = {), (7) 
pdpV {s -ip)v - {m^\)w = (8) 

We eliminate w and, after introducing a new independent 
variable z = —2ikp^ obtain a hypergeometric equation 



zv" + (2s + 1 - z)v' -{s-ip)v = 0. 
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The solution regular at 2; = is given by the confluent 
hypergeometric function [1^ : 



v{z) = A iFi {s - ip, 2s + 1, ^) , 
s — ip 



(10) 



w{z) = A ^ iFi (s + 1 - iP, 2s + 1, ^) , (11) 

where A is a normalization factor. The expression for w 
was obtained using Eq.([8|) and an identity for iFi [14 1. 



The solution ([5]) , ([TO]) , ([II]) of the Dirac-Kepler problem, 
regular at p = 0, can be used to evaluate the polarization 
charge in the subcritical case \p\ < 1/2. This can be 
done most easily using the scattering phases, given by 
the behavior of w and v at large p. Using the asymptotic 
form of the functions iFi [l5|, we find 



^^if3 ln(2/cp) A* e"*'^ ln(2/ep) ^-2ikp 

v{p)= ,o/..As ^ ^ip) = 77Tr:x; ' (12) 
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where the parameter A depends on m and P but not 
on k. From (p!2|) we see that p^~^w{p)ex.p{ikp) and 
p^~^v{p) ex.p{ikp) indeed describe the incoming and the 
outgoing waves, characterized by relative phase 

v/w = e2^^-(^)+2^^^, Sm{k) = p\n{2kp) + arg A. (13) 



The log dependence in Eq. (p!3|) is typical for the phase 
coming from 1/r Coulomb tail [12]. 

The scattering phases Sm{k) can be used to find 
the polarization charge pulled on the origin. How- 
ever, a straightforward application of the Friedel sum 
rule |l3 , involving phases evaluated at the Fermi 
level, encounters a difficulty due to the position and en- 
ergy dependence of the phases in Eq. (p!3|) . Since this 
may indicate that the polarization is distributed rather 
than localized at p = 0, we proceed with caution. 

To evaluate the excess particle number Qpoi(p) in the 
interval < < p, we note that the states with wave- 
length greater than p, i.e. \k\ < 1/p, contribute negligibly 
to Qpo\{p). Thus we can write the sum rule [Ifi, J/7.] as 



Qpoi(p) 
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where the minus sign corresponds to that in Eq.(|6]). Con- 
veniently, the expressions (p!3|) for Sm{k), valid at kp > 1, 
can be used to evaluate ([E]). However, since Sm{k) de- 
pend on the product kp, they yield a p-independent result 
for Qpo\{p)- We therefore conclude that the polarization 
charge is concentrated on the lattice scale p ^ tq. 

To independently verify the conclusion about polariza- 
tion at p < 1/2 concentrated at p < ro we evaluated it 
directly using the eigenstates given by Eqs. ([5|) . (p!Q|) . (pTj) . 
This calculation involves an energy cutoff introduced at 
the bottom of graphene band, corresponding to k ^ 
. We found nonvanishing contribution to polarization 
charge only on the cutoff scale, leading to an expression 
^poi(p) = —Qp^ip) (the second term in Eq.([2])). This form 
of polarization charge can be independently justified by 
the RPA method giving qp = ^p. Our numerical 

results at < < 1, presented in FigiH yield a very 
similar dependence, nearly linear at \P\ < ^, indepen- 
dently confirming the above analysis. 

The behavior of the scattering phase and of the po- 
larization charge changes when the potential strength \p\ 
exceeds |m + ^| for one or several values of m. For such 
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FIG. 1: Local polarization charge found from numerical solu- 
tion of a tight-binding problem on a honeycomb lattice. The 
charge was placed in the middle of a rectangular region of size 
ni X 712 at the lattice plaquette center or edge center (see leg- 
end) . Shown is the net polarization at a distance of less than 
5 lattice constants from point charge, which agrees with the 
prediction of a perturbative RPA calculation (dashed line) . 



supercritical P Eq. (p!Q|) is not the only possible solution. 
Adding another solution of the equation ([9]), we write the 
function v{z) in the form 

v{z) = A iFi (i(7 -/?),! + 2^7, z) (15) 
+ 5 z-2^^iFi(-i(7 + /?),!- 2^7, z), 



where 7 = Ims 



and z 



-2ikp. 



With the help of the relation [14] we find 



w{z) =-iAr]iFi{l^ij-ip,1^2ij,z) (16) 
-i{B/r])z-^'\Fi (1 - i7 - iP. 1 - 2^7, 



where r] = y fqp^- 

To find the relation between A and B we consider our 
solution at small distances, p ^ ro <C l/k: 



v{p) 
w{p) 



-iAf] — i{B/r])e 



-2i7 ln(2/cp) 



(17) 
(18) 



Without loss of generality we use the boundary condition 
'02(p = ^0) = O5 which enforces zero wavefunction on one 
of the graphene lattice sites. (Similar boundary condition 
was used to describe graphene zigzag edge in the Dirac 
equation framework [18].) Solving the equation v{ro) = 
w{ro)^ we find the relation 
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The product kro is very small for typical /c, making the 
phase factor e*^^^^ a rapidly oscillating function of k. 

To better understand the role of the second solution, 
let us take another look at the subcritical case, < |m+ 
^1, when the parameter 5, Eq.(|6]), is real. In this case, two 
independent solutions are still provided by Eqs. (p!5|) . (p!6|) . 
whereby is replaced by s. Applying the boundary 
conditions the same way as above, instead of (p!9|) we find 



B/A(x (2/cro )^* <^ 1, which indicates that the second 
solution plays no role in the subcritical case. 

To link x(^)? Eq. (p!9|) . to the scattering phase, we write 
our solution for v{p)^ ^{p) at large distances pk 1, 
again using the asymptotic expression for iFi [15|], 
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where ^^3,7 = T(l^2i-f) /T{l^i-f^i(3). We note that (f2Q|) 
automatically satisfies the current conservation require- 
ment \v\ — \w\. The relative phase of v and defined 
as in Eq.([T3j), thus equals 



^m(^) = ^(^) +/31n2/cp + arg^/3,^ 



(21) 



where 



d(k) — arg 
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The last two terms of the phase ([2T]) are identical in 
form to (p!3|) . They represent spherical wave "deformed" 
by the tail of Coulomb potential at large distances and, 
as we discussed above, give no contribution to polar- 
ization charge at finite p. The term 0(k\ however, 
arising from the boundary condition at small p via the 
phase X — arg(5/A), Eq.([T9|), makes the behavior com- 
pletely different from that found for < |m + ||. 

The phase d(k) dependence on /c, arising via the phase 
x(A:), is quite peculiar. To understand the relation Q vs. 
X, we find the winding number of 0{k) that depends on 



a = 
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< 1 if /3 > 0, 
> 1 if /3 < 0. 
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(We recall that < 7 < \/3\.) The phase winding is 
thus controled by the first term of ([22|) at /3 > and by 
the second term at /3 < 0, allowing us to write it as 

e{k) = - sign/3^ + A(9(/c), (24) 

where is an oscillatory periodic function of x- 

It is instructive to compare the behavior of AO{k) at 
strongly overcritical and nearly critical /3. For large \/3\ 
we have 7 and |a| ~ e~^^^, and thus the oscillatory 
part is exponentially small: 

AO{k) ^ sign/?e-^^^ sin(arga - x{k)). (25) 

In the opposite limit of nearly critical P ~ pc = ^l^ + ^l^ 
we expand \a\ in 7, which is small in this case, to find 



|a| = l-e, ? 



27r7 



+ 0(7')- (26) 
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In this case 6*, as a function of x, is a staircase 

6* = ^ - arctan U tan ^^y^ ] , C = arg a, (27) 



with steps of height tt, width 27r, and corners rounded 
on the scale ^. The staircase slope — ^ sign/3 corresponds 
to the first term in Eq. ([M|) . The oscillatory part AO{k) 
manifests itself in the local density of states around the 
impurity (see peak in the e < region in Fig l2] inset). 

To analyze the contribution of the phase 0{k) to the po- 
larization density, we suppress the periodic part A^. Us- 
ing the expression ([22]) we find 0{k) = — 7sign/31n2/cro. 
Substituting it in the Friedel sum rule, Eq.([T4j), we find 

qMp) = -ivfclM = _ ,ign/jl^ m JL. (28) 

TT TT 2ro 

From npoi(p) = {27Tp)~^dQpoi/dp^ we find the polariza- 
tion density ([2|). When the parameter s is complex in 
more than one channel, one has to consider a sum 

, , TV sign /3 r~ ~. 772 

^poi(p) = E yf^ - + ^) • (29) 

|m+i|<|/3| 

For large P 1, replacing the sum by an integral, we 
recover the expression npoi(p) = A^/3|/3|/(47rp^) found in 
0] by the Thomas-Fermi method. 

For supercritical /3 > ^, as discussed above, the scat- 
tering phase becomes sensitive to the physics at small 
distances, p ^ ro. This leads to pronounced interference 
of the incoming and outgoing waves, which is manifest in 
the local density of states (LDOS), 

K^'^) = J^EW^^'^)^ ^^ = -)^' (30) 

m 

with an appropriate normalization of the two-component 
wave function ip^ Eq.([5]). We evaluated the sum in ([30]) 
numerically, using the expressions (p!5|) and (p!6|) . For 
< \(3\ < ^ LDOS does not deviate too much from 
1^/3=0 kl (see Fig 12] inset). For supercritical /3, however, 
LDOS develops pronounced oscillations in both position 
p and energy k^. The crossover from a non-oscillatory to 
oscillatory behavior dit p ^ ^ becomes sharp at p ^ tq. 

The standing waves in LDOS (|3Q]) at /3 > ^ are dif- 
ferent from Friedel oscillations, since kp = in our case 
(Fermi level at the Dirac point). As illustrated in Fig 121 
their spatial period scales inversely with energy, so that 
the maxima occur at k^p = (n + ^)7r, which is simi- 
lar to the oscillations in LDOS studied in carbon nan- 
otubes [l^ . As in Ref. [l^ , the energy-dependent spatial 
period can be used to obtain direct information about 
Fermi velocity in graphene. 

The spatial structure predicted around supercritical 
Coulomb impurities, which extends up to the vacuum 
polarization cloud size p^, will be affected by finite tem- 
perature T > = hvp/pi,^ and also by carrier doping 
away from neutrality by Sn > strong enough to in- 
duce screening at distances less than p^. 

To summarize, we found that the excess charge — ^ of 
supercritical impurities in graphene is fully screened by 
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FIG. 2: Standing wave oscillations in LDOS, Ea.(l30l). vs. en- 
ergy scaled by the distance from impurity p. Maxima (min- 
ima) occur at half-integer (integer) values of kp/n. LDOS is 
shown for an over critical /3 = 0.6 and several values of given 
in the units of 10"^ ro- Inset: Oscillations in LDOS, appearing 
for 1/3 1 > I on the top of the unperturbed density of states, 
which is subtracted in the main figure (sp — hvp/ p)- 



the Dirac vacuum polarization. The large screening cloud 
size and the standing wave oscillations predicted within it 
can be directly probed by STM technique. The sharp de- 
parture fom linear screening for supercritical impurities 
represents an interesting example of nonlinear screening 
that can be realized in graphene. Our estimates for the 
critical charge, using /^rpa ^ 5, yield an experimentally 
convenient value Zc ~ 1, making experimental tests of 
these effects in graphene practical. 
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